library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # TS padding
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz
library(forecast) # TS library
library(TTR) # SMA function
library(tseries) # adf.test
# supaya semua plot memiliki theme_minimal()
theme_set(theme_minimal())Gunakan data souvenir untuk melakukan fitting model dan forecasting menggunakan metode exponential smoothing. Data tersebut merupakan data penjualan souvenir per bulan selama 7 tahun dari bulan Januari 1987 sampai Desember 1993.
Load data menggunakan fungsi scan() untuk emmbaca data menjadi sebuah vector, lalu lakukan visualisasi dengan base plot()
# import data
souvenir <- scan("data_input/fancy.dat")
# line plot
plot(souvenir, type = "l")Buatlah object Time Series ts(data, start, frequency), lalu visualisasikan dengan autoplot()
# object TS
souvenir_ts <- ts(data = souvenir,
start = c(1987, 1),
frequency = 12)
# autoplot
souvenir_ts %>%
autoplot()Insight: Apakah TS termasuk tipe additive/multiplicative? multiplicative
Lakukan decomposition sesuai dengan tipe additive/multiplicativenya, lalu buatlah bar plot untuk menganalisa komponen seasonalitynya.
# decomposition
souvenir_dc <- souvenir_ts %>%
decompose(type = "multiplicative")
souvenir_dc %>%
autoplot()# seasonality analysis
data.frame(
Month = factor(month.abb, levels=month.abb),
Seasonality = souvenir_dc$seasonal
) %>%
distinct() %>%
ggplot(aes(x = Month, y = Seasonality)) +
geom_col()Insight: Bagaimana trend dan seasonality dari penjualan souvenir?
Lakukan train-test splitting, dimana 1 tahun untuk data test dan sisanya untuk data train
# data test
souvenir_test <- tail(souvenir_ts, 12)
# data train
souvenir_train <- head(souvenir_ts, -12)Buatlah minimal 2 model exponential smoothing dengan fungsi berikut:
HoltWinters()ets()Note, pemilihan model berdasarkan:
souvenir_hw <- HoltWinters(souvenir_train, seasonal = "multiplicative")
souvenir_hw_add <- HoltWinters(souvenir_train %>% log(), seasonal = "additive")
souvenir_ets <- ets(souvenir_train, model = "ZZZ")Gunakan forecast() untuk masing-masing model yang Anda buat sebelumnya:
souvenir_hw_f <- forecast(souvenir_hw, h = 12)
souvenir_hw_add_f <- forecast(souvenir_hw_add, h = 12)
souvenir_ets_f <- forecast(souvenir_ets, h = 12)Visualisasikan data actual (train dan test) dengan hasil forecasting
p <- souvenir_train %>%
autoplot(series = "Train") +
autolayer(souvenir_test, series = "Test") +
autolayer(souvenir_hw_f$mean, series = "HW Multiplicative") +
autolayer(souvenir_hw_add_f$mean %>% exp(), series = "HW Additive") +
autolayer(souvenir_ets_f$mean, series = "ETS")
plotly::ggplotly(p)accuracy(souvenir_hw_f, souvenir_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 189.1763 2200.411 1512.406 -0.1299785 13.13863 0.4038682
#> Test set -4060.1992 4541.688 4060.199 -21.1511565 21.15116 1.0842231
#> ACF1 Theil's U
#> Training set 0.1829319 NA
#> Test set 0.5148249 0.6278038
accuracy(souvenir_ets_f, souvenir_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 84.47865 1805.020 1228.946 -0.5220307 11.80610 0.3281741
#> Test set -1019.41169 3790.485 3389.888 -10.2903652 18.22484 0.9052254
#> ACF1 Theil's U
#> Training set 0.2111953 NA
#> Test set 0.7158307 0.5309041
# cara mengevaluasi apabila data multiplicative diperlakukan sebagai additive
accuracy(souvenir_hw_add_f$fitted %>% exp(), souvenir_train) # training#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -55.25421 2302.798 1552.995 -1.715425 13.25152 0.1180214 0.3660692
accuracy(souvenir_hw_add_f$mean %>% exp(), souvenir_test) # testing#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -5769.926 6138.801 5769.926 -25.9165 25.9165 -0.04706218 0.6797808
Kesimpulan: Berdasarkan nilai error di atas, model yang terbaik berdasarkan metric MAPE adalah souvenir_ets yaitu ETS(M,A,M).
Bagian ini opsional, namun penting untuk diketahui
Misalkan selanjutnya ada kebutuhan untuk melakukan forecasting periode Jan 1994 - Dec 1994 (satu tahun kedepan dari data souvenir). Maka practice yang baik, biasa kita re-train semua data kita menggunakan best model pada tahap sebelumnya
# misal pada tahap sebelumnya, model terbaik adalah dengan fungsi ets() model ZZZ
souvenir_ets_full <- ets(souvenir_ts, model = "ZZZ")
souvenir_ets_full_f <- forecast(souvenir_ets_full, h = 12)
# visualisasi
souvenir_ts %>%
autoplot(series = "Full Data") +
autolayer(souvenir_ets_full_f, series = "Forecast 1 year")Model souvenir_ets hanya di train dari Januari 1987 sampai Desember 1992, sedangkan souvenir_ets_full di train dengan keseluruhan data (sampai Desember 1993).